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It is known chat physical conditions in the solar atmosphere are 
strongly controlled by solar magnetic fields. Therefore, the measurement 
of solar magnetic fields is essential to the understanding of the physics 
of the solar atmosphere . The ground-based measurements have been 
obtained for some time. Making solar magnetic field measurements in 
space would have better resolution. In this study we conducted an 
investigation and assessment and planning to build such a facility on 
board space station (or a moon-based observatory) to measure solar 
magnetic fields from space. This was accomplished through participation 
in scientific studies. Specifically the following were completed; 


Dr. Rainer Kress, visited MSFC/Space Science Laboratory and UAH 
to discuss and gave seminar on mathematical methods related to 
extrapolation of solar magnetic fields March 15 - 20, 1991. 

These discussions led to an article "A Comparison Between 
Progressive Extension Method (PEM) and Iteratiave Method (IM) for 
Magnetic Field Extrapolations in the Solar Atmoshere" by S. T. 
wu, M. t. Sun and T. Sakurai, published in Mem S.A.It., 1990, 

Vol. 61, No. 2, pp 477-484. (See attachment) 

Dr. Petrus Martens visit UAH and MSFC/Space Science Laboratory on 
November 27-30 1990 to discuss the shear motion instability 
related to flare onset. These discussions lead to an article to 
be published entitled "Shear- Induced Instability and Arch 
Filament Eruption; A Magnetohydrodynamic (MHD) Numerical 
Simulation" by S . T. Wu, M. T. Song', P. C. H. Martens, and M 
Dryer, Solar Physics, 134, 353-377, 1991. (See attachment) 

Dr. J. P. Rozelot visited UAH and MSFC/Space Science Laboratory 
on November 26 - 28, 1990 to discuss the reconstruction of the 
solar cycle and make brief presentation on large active mirrors 
in aluminum, for new generation of telescopes. 


These discussions and seminars provided Che advancement of kev 
concepts and technology in this planning for observations of solar 
magnetic fields with space Qualified instrumentation. 
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A COMPARISON BETWEEN PROGRESSIVE EXTENSION METHOD (PEM) AND 
ITERATIVE METHOD (IM) FOR MAGNETIC FIELD EXTRAPOLATIONS 
IN THE SOLAR ATMOSPHERE 


S. T. Wu and M. T. Sun 

Center for Space Plasma and Aeronoraic Research 
and Department of Mechanical Engineering 
The University of Alabama in Huntsville 
Huntsville, AL 35899 USA 
and 

Takashi Sakurai 

National Astronomical Observatory 
Mataka, Tokyo, 181 Japan 


ABSTRACT 

In this paper we present a comparison between two numerical 
methods for the extrapolation of nonlinear force-free magnetic 
fields viz. (i) the Iterative Method (IM) and (ii) the 
Progressive Extension Method (PEM) . The advantages and 
disadvantages of these two methods are summarized and the 
accuracy and numerical instability are discussed. On the basis 
of this investigation, we claim that the two methods do resemble 
each other qualitatively. 


I. INTRODUCTION 

It is wellknown that the magnetic fields play a dominant 
role in all physical features which appear in the solar 
atmosphere; for example, the observed filamentary structures in 
the chromosphere seen in H^ (Martin, 1980), and coronal loops 
seen in UV (Cheng, etal. 1982) and X-rays (Antonucci etal. 1982; 
de Jager et al. 1983). All these structures in the solar 
atmosphere are generally considered to be aligned along the 
magnetic field (Zirin, 1971; Poletto, et al., 1975). 
Physically, these structures can be interpreted as plasma 
confined by the magnetic field- Hence, a detailed and 
quantitative analysis of these structures require a 
quantitative knowledge of the magnetic field in the solar 
atmosphere. Presently, measurements of magnetic fields are 
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confined to the pho.jspheric level; therefore, _n higher levels 
(i.e. chromosphere and corona) the magnetic field can only be 
obtained through numerical extrapolation using the measured 
photospheric magnetic field as the source surface, as 
demonstrated in the early work of Schmidt (1964) , Altschuler and 
Newkirk (1969), Nakagawa and Raadu (1972). All these early 
extrapolation methods are restricted to the linear 
approximation, which physically represents current~free field 
(potential field) or constant current-to-magnetic field ratio 
(linear force free field) . It has been shown that these 
representations are far from realistic in describing the 
observed features in the solar atmosphere (Schmahl et al., 
1982) . 

In order to improve our understanding of the physical 
structures of the solar atmosphere it is necessary to have 
quantitative knowledge of the magnetic field. Therefore, a 
number of extraplation methods is developed to meet the demands. 
The mathematical model using a force free configuration on the 
basis for the extrapolation of photospheric vector magnetograms 
to obtain the coronal field has been given by Aly (1989) and Gary 
(1990). In particular, Gary (1990) presented an excellent 
summary and assessment on the present available extrapolation 
methods from a theoretical point of view. In this paper, a 
comparison between the progressive extension method (PEM) and 
iterative method (IM) is presented. The rationale for choosing 
these two extrapolation techniques for comparison is that they 
are based on observed photospheric level fields and have 
practical applications. A brief description of the theoretical 
background of these two techniques is presented in Section 2 . 
Numerical results of direct comparison are included in Section 
3 . The discussion of advantages and disadvantages of these two 
techniques and their possible physical consequences are 
presented in section 4 . 


II. THEORY AND TECHNIQUES 

On the assumption of magnetohydrostatic equilibrium in the 
solar atmosphere, the mathematical model describing such an 
equilibrium state may be written as 

-9p +JxB-pg=0, (i) 

where p is the hydrostatic pressure and will be represented by 
the equation of state, 

p = pRT , (2) 

with p and T being the mass density and temperati^re respectively . 
The other symbols have their usual meanings; B is _the magnetic 
field and J, the current density, is related to B by 
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J 


7 X B . 


(3) 


Finally, g is the gravitational acceleration. Physically, 
there are three different orders of approximation to determine 
the magnetic field configuration. The first and second order 
approximations are the current free (potential) and force-free 
magnetic field, respectively. Within these orders of 
approximation the magnetic force vanishes, and the pressure 
force is balanced by the gravitational force which leads to the 
hydrostatic equilibrium in the solar atmosphere. Under these 
circumstances, the mathematical model for the magnetic field 
configuration can be represented by 

7 X B = aB , (4) 

This expression possesses three different physical meanings, 
which are: (i) a = 0, corresponds to the current free case in 
which the magnetic field is potential, (ii) a = constant, 
corresponds to the linear force-free magnetic field which 
implies a constant current-to-magnetic field ratio in a region 
and (iii) a = a(r) , corresponds to the nonlinear force-free field 
which implies a non-constant current-to-magnetic field ratio in 
a region. 

Finally, the third order of appoximation is the 
magnetohydrostatic equilibrium in the solar atmosphere which is 
given by Eq. (1). If there is information on_B and p on the 
source surface, it is possible to extrapolate B and p upward. 
Since there only are measurements of the magnetic field on the 
source surface (photosphere) , it is not possible to extrapolate 
magnetohydrostatic equilibrium field-configurations at the 
present time. 

In the meantime, we shall focus our attention on the 
nonlinear force-free field configuration. For the purpose of 
this paper, we have selected two techniques for this 
investigation. These tv;o techniques are progressive extension 
method (PEM) (Wu et al., 1985, 1990) and iterative method (IM) 
(Sakurai, 1981) . A brief description of these two methods is 
presented below: 

Progressive Extension Method (PEM) 

The progressive extension method is formulated as an 
initial-value problem (i.e., Cauchy problem) using a finite 
difference scheme which is similar to a Taylor expansion. A 
detailed description of this method is given by Wu et al . ( 1990) . 
They have demonstrated the usefulness of this method, and the 
numerical algorithm has been verified by extrapolation of an 
analytical solution (Low, 1982) . 
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Iterative Method ( 


A number of authors (see references in Gary, 1990) have 
utilized an iterative method originated by Grad and Rubin (1958) 
to extrapolate the nonlinear force-free magnetic field from 
boundary data. For convenience, we simply choose the iterative 
method developed by Sakurai (1981) in this study. His method is 
based the integral equation representation of Eg. (1), and the 
discretization is made by the technique of finite element 
method. A detailed description of this technique was given by 
Sakurai (1981), and we shall not repeat it here. 


III. NUMERICAL RESULTS 

In order to make comparison between the PEM (Progressive 
Extension Method) of Wu et al. (1985, 1990) and the IM ( Iterative 
Method) of Sakurai (1981) , we have chosen the vectoral magnetic 
field observed at Okayama Astrophysical Observatory on May 26, 
1985 (Sakurai and Makita, 1986) as the boundary for 
extrapolation using these two methods. The observed magnetic 
field vector is shown in Figure 1. 


MC526C DATE: 83/5/26 TIME (JST): 10 2 13 -11 13 5 
OBSERVED FIELD VECTOR 



Figure 1. Magnetic field vector observed at Okayama 
Astrophysical Observatory on May 26, 1983. Solid 
and dotted contours show positive and negative 
longitudinal fields, respectively, with levels ± 10, 
20, 50, 100, 200, 500 G. Arrows indicate the 

transverse vector. 
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Usinq these obst. '/ational data as a sourct surface, we 
obtained the nonlinear force-free field configuration by using 
the above mentioned two methods as shown in Figure 2^ where 
Figure 2a is obtained by using the IM and Figure 2b by using PEM. 
In addition we have extrapolated the potential 
configuration using PEM in comparison with the potential field 
given by Sakurai and Makita (1986) , see Figure 3. From these 
results, observe that the deduced magnetic field configurations 
albeit not identical, in fact, qualitatively resemble each other 
to a large extent. 



Figure 2. Nonlinear force-free field lines computed by (a) 
Iterative Method (IM) and (b) by Progressive 
Extension Method (PEM) using the data shown in Figure 


1 . 



Figure 3. (a) Potential field lines computed by IM and (b) 

potential field lines computed by PEM using the 
observation given in Figure 1. 
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IV. DISCUSSION 


Before we analyze the causes of these differences seen in 
the two extrapolations we review the fundamental differences 
between the two methods. These differences can be summarized as 
follows : 

1. The Iterataive Method (IM) specifies the value of a on a 
portion of the boundary plane (e.g. on a positive field 
region) and cannot assign the value of a on the whole 
boundary plane, since that would introduce an 
inconsistency in the extrapolation process. The 
values of a in the whole boundary plane are determined 
by the observed data for PEM. In this fashion, there 
is an electric current only along the particular field 
line in the IM extrapolation, while the electric 
current is distributed in the whole domain of 
calculation for the PEM extrapolation. 

2. The IM type of extraplation is convergent only for 
small values of a. Physically, this implies that the 
electric current in the region of interest must be 
small. On the other hand, the PEM type of 
extrapolation does not have this limitation. 
However, the accuracy of the computed a-value 
deteriorates at the points near the neutral line (i.e. 
Bj ^ 0) . This may cause a misrepresentation of the 
magnetic field configuration. The grid size of the 
extrapolation is controlled by the numerical stability 
criteria as given by Wu et al. (1990). 

3. The fact that the value of a is assigned at one of the 
two foot points of a particular field line in the IM 
while the values of a are determined on the entire 
boundary surface in the PEM makes it difficult to match 
and compare the field lines for these two different 
methods . 

On the basis of these differences of extrapolation 
procedures, we may understand why the magnetic field 
configurations obtained from the same data with these two 
methods are not identical. For example. Figure 2, shows some 
differences in magnetic field-line configurations, but the 
lines connecting different regions of polarities are quite 
similar. Note that for two regions of opposite polarities near 
the right center, the PEM extrapolation doesn't show any 
connection by field lines, while the IM type extrapolation does. 
However this is due simply to the fact that the field lines in 
this region are very low and short, and cannot be discerned in 
this drawing. Plots of the front view of Figure 2b, clearly 
indicate that the regions are connected by field lines (marked by 
A) as shown in Figure 4 . 
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Figure 4. The front view of the nonlinear force-free field 
computed by PEM using the observtion given in Figure 
1. It should be noted that the field lines near the 
top are not accurate due to numerical procedure as 
discussed by Wu et al. (1990). 

We further notice that the configuration of the field lines 
obtained by IM extrapolation is very similar to a potential field 
line configuration. This is because the IM requires that the 
value of a be small (i.e. slightly deviating from potential) . 
On the other hand, the PEM extrapolation does not have this 
limitation. It is understood that the degree of deviation from 
a potential field depends on the value of a, that is the strength 
of the local electric current. Therefore, the configuration of 
magnetic field lines is affected. 

In summary, we conclude: 

(i) Both methods do produce qualitatively similar results, 

(ii) The accuracy of PEM has been verified by an analytical 
solution (Wu et al. 1990) ; verification of IM is still 
needed. 

(iii) There are limitations on the value of a for IM, but not 
for PEM. 

(iv) The accuracy for PEM deteriorates when the height of 
extrapolation exceeds one third the horizontal length, 
because of the propagation of the accummulated 
numerical errors at each level (Wu et al. 1990). 
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SHEAR-INDUCED INSTABILITY AND ARCH FILAMENT 
ERUPTION: A MAGNETOHYDRODYNAMIC (MHD) 
NUMERICAL SIMULATION 
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Abstract. We investigate, via a two-dimensional (nonpianar) MHD simulation, a situation wherein a bipolar 
magnetic field embedded in a stratified solar atmosphere (i.e., arch-filament-like structure) undergoes 
symmetrical shear motion at the footpoints. It was found that the vertical plasma flow velocities grow 
exponentially leading to a new type of global MHDMnstability ih3X could be characterized as a ‘Dynamic 
Shearing Instability’, with a growth rate of about 8 V^a, where is the average Aifven speed and a ~ ’ 
is the characteristic length scale. The growth rate grows almost linearly until it reaches the same order ot 
magnitude as the Aifven speed. Then a nonlinear MHD instability occurs beyond this point. This simulation 
indicates the following physical consequences: the central loops are pinched by opposing Lorentz torces, 
and the outer closed loops stretch upward with the vertically-rising mass flow. This instability may apply 
to arch filament eruptions (AFE) and coronal mass ejections (CMEs). 

To illustrate the nonlinear dynamical shearing instability, a numerical example is given tor three different 
values of the plasma beta that span several orders of magnitude. The numerical results were analyzed using 
a linearized asymptotic approach in which an analytical approximate solution for velocity growth is 
presented. Finally, this theoretical model is applied to describe the arch filament eruption as well as CMEs. 


1. Introduction 

More than a quarter century ago. Gold and Hoyle (1960) suggested that horizontal 
photospheric motion can move the footpoints of magnetic field lines and twist the flux 
tubes because of the highly electric conducting plasma at the photospheric levels. A 
number of investigators (Tanaka and Nakagawa, 1973; Low and Nakagawa. 1975; 
Low, 1977; Klimchuk, Sturrock, and Yang, 1988: Klimchuk and Sturrock, 1989) 
studied the evolution of force-free fields and its role in energy storage (build-up) for solar 
flares. 

* Permanent address: Purple Mountain Observatory, Nanjing, China. 
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All of these studies were limited to the case of magnetostatics: self-consistent dynami- 
cal effects were ignored. Recently. Wu, Hu. and Nakagawa ( 1983). Wu. Hu. and Krall 
(1984). and Wu ei al. (1986) presented a self-consistent MHD model for the purpose 
of e.xamining flare energy build-up and wave-mass interactions due to shear and 
converging-diverging motions at the photospheric level. More recently. Mikic. Barnes, 
and Schnack (1988) and Biskamp and Welter (1989) have presented numerical results 
on the dynamical evolution of a magnetic arcade type due to shear motion. However, 
their models are restricted to symmetric boundary conditions, while in this stud\ 
self-consistent boundary conditions were used (see. for example. Wu and Wang, 1987 ; 
Nakagawa. Hu. and Wu. 1987). 

In this paper, we use the time-dependent MHD simulation model devised by Wu. Hu. 
and Nakagawa ( 1983) to reveal a nonlinear solution for the evolution of the magnetic 
field configuration driven by shear motion. In this solution, we find that the plasma 
velocity in the vertical plane perpendicular to the shear, grows exponentially in a process 
which can be analytically described by a linear MHD instability. This upward velocity 
steadily increases until it reaches the average Alfven speed. .\t later times, a nonlinear 
instability sets in. .A field line pinch occurs in the lower shear region in the numerical 
results. .At the same time, mass and field line expulsion appears in higher parts of the 
region and the closed field tends to open locally. We suggest that these new effects (i.c.. 
mushroom cloud-like flow, pinch, and expulsion) can explain the formation of current 
sheets, the opening of a closed bipolar field, and the ability of particle streams to escape 
from the solar surface. Specifically, we suggest that this model applies to the eruption 
of arch filament systems (AFEs) and their relation to non-flare-associated coronal mass 
ejections (CMEs). The mathematical description of the model and numerical results arc 
given in Section 2. .A general physical interpretation of these results is presented in 
Section 3. .An application of this model to specific coronal phenomena is given in 
Section 4. and the concluding remarks are presented in Section 5. 

2. Numerical Simulation 

In order to illustate how shear-induced non-equilibrium occurs, we use a theoretical 
model in which a two-dimensional bipolar field undergoes a steady shear velocity at the 
footpoints of its magnetic loops. The shearing motion is sketched in Figure 1(a), and 
the initial bipolar field is shown explicitly in Figure 1(b). 

First, we perform a simulation of the dynamic response of the bipolar field to the 
shear. Then we use an analytical method to interpret the simulation results. The 
simulation model is based on a two-dimensional, time-dependent, MHD model (Wu. 
Hu, and Nakagawa, 1983: Hu and Wu, 1984) with an improved FICE (Full-Implicit- 
Continuous-Eulerian) numerical scheme (Wu and Wang, 1987). Symmetrical side 
boundary conditions have been replaced with non-reflecting boundary conditions. 
This implies that the physical phenomena are determined by the solution at a specific 
time and are not determined by the specified boundary conditions as in the case 
studied by Mikic. Barnes, and Schnack (1988). The physical conditions on these 
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(b) 

Fig, 1. (a) Sketch of a two-dimensional bipolar magnetic field that is subjected to a footpoint shearing 

motion as indicated by the arrows, (b) Explicit bipolar magnetic topology prior to the shearing motion (see 
Equation ( 1 )). The photospheric boundary extends to ixt = 8.4 x 10^ km in both directions from the origin. 
The vertical extent into the corona is to y = 8 x 10^ km. The positions y = Vj. Vj, — ,Vjn indicate the 
vertical levels at which horizontal surveys will be shown of various physical quantities during the shearing 

motion at the footpoints. 


two side boundaries are determined mathematically through compatibility relations 
that are given in detail by Wu and Wang (1987). Thus, the computation domain 
(i.e., I.xl <8.4 X 10^ km, 0<y<8 x 10^ km) consists of three free non-reflecting 
boundaries (i.e., top and sides), while the bottom boundary (>’ = 0) is treated with the 
method of projected characteristics (Nakagawa, Hu, and Wu, 1987 ; Hu and Wu, 1984). 
The basic equations for this model are the time-dependent MHD equations with infinite 
conductivity, no viscosity and symmetry in one direction (Wu, Hu, and Nakagawa. 
1983). Solar gravity, plasma pressure gradients, and compressibility are explicitly con- 
sidered. None of these characteristics were considered in the work of Mikic, Barnes, 
and Schnack (1988), and Biskamp and Welter (1989) have only considered compressi- 
bility in a special way. 
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The initial conditions are (see Figure 1(b)): 


A exp( - ■ )- 

^ rtJ 

T, = T. , ty, ty, ly = 0 . 



= Br, [cos(a.v)] e 

. 5, 0 = - 5o [sinia.x)] e? ' " . 

5,0 = 4 . 

( 1) 


a = 7r/2.V|, . ,v,| = 8.4 x 10^ km . g = 2.71 x lO'* cm s ^ . 

The plasma parameters are taken to be = 1.67 x 10“ g cm ‘ ^ and = 10^ K. The 
scale height ' = RT^jg ^ 6.1 x 10^km)anda“‘ ~ 6.3 x 10^ km are the same order 
of magnitude. These parameters are representative for solar conditions at the higher 
chromosphere and lower corona. The computation grid points are: 

.V, = - 8.4 X 10^ + (/ - l)d,x . /■ = 1, 2 22 . 

>v = (7-l)4v. 7 = 1.2 11 , 

zl.x = Ay = 8 X 10" km ~ 1 arc sec . 

The non-reflecting boundary conditions, as noted above, are used for the top ( r = r, , ). 
left-hand side (.x = .v , ), and right-hand side (.x = The conditions at the bottom 
boundary ( r = r,) are taken as follows: 

P = Pc- , T = T^, B, = 5,.o , = 0 . but ty, v. ^ 0 . 

sln(a.x) if kxi < 5.2 x 10^ km . 

(6.8 X 10" - .Vi) 

IV. (sgn.x) sin(5.2 X 10"a) 

1.6x10" ^ (2) 

if 5.2 X 10" < i.xi < 6.8 X 10" km . 

0 if 6.8 X 10" < ..xj <8 X 10" km . 

The other physical quantities (p. T. l\„ 5.) are computed by means of the 

compatibility equations for the non-reflecting boundary condition which assures the 
consistency of the numerical computation. 

In order to understand the general physical behaviour of the nonlinear solution from 
the mathematical model, we have performed three numerical experiments. These three 
cases use combinations of magnetic field intensity and magnitudes of the shear velocity. 
The results for these three cases are described as follows. 

2.1. Large plasma beta (jSg ~ 154) 

In this numerical experiment, we choose the initial plasma beta ()?o) to be 154 where 
is defined as = Pn/(^o/8^) with p„ and B^, being the plasma pressure and magnetic 
field strength at the lower boundary (i.e., v =Vi). This is not a physically realistic 
case for a solar active region; but it does provide a basis for comparison with the other 
cases. This case corresponds to a local, exceedingly low, magnetic field strength of 
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2.12 G at the origin, .x = r = 0. as shown in Figure Kb). The shear velocity, a; , was 
taken to be 5 km s ' ' . Figure 2 shows the evolution of the magnetic field lines due 
to the shear motion at 200 s < r < 3200 s. It is useful to examine the ev'olutionarv 
behaviour at various Alfven limes (defined as = [Jv (or d.v)],' F^ ~ 1700 s where 


200s 



600s 





2000s 


2400s 


(0 


2800s 


(g) 


3200s 






p. ~ 15 4, 1700s 

Fig. 2. Magnetic field line evolution as a function of time during induced footpoint shearing motion for 
case (i). Pq = 154 and the Alfven time, = 1700 s. The horizontal axis represents the distance from 

x^ ... X 22 as shown in Figure 1(b). 
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^ \ = ^n!\ -^^Po = km s ‘ ' ). During the early stages of evolution (that is. within the 
first Alfven time), the magnetic field lines rise together in an orderly fashion in response 
to the shearing motion. This behaviour is also presented in the analytical solution of Low 
( 1 98 1 ) and the force-free numencal solutions of Klimchuk and Sturrock ( 1 989) although 
they do not consider dynamics and gravitational effects. .After the first .Alfven time 
period, the evolutionary behaviour of the field lines becomes more complicated. 
Nonlinear interactions take place between the shear-induced mass motion, maanetic 
field and gravity with the result that in some regions the field lines are bunched loaether 
to form a current sheet (see Figures 2(g) and 2(h)). Further understanding of these 
phenomena is provided by the representation of the shear induced mass motion as 
shown by the vectorial velocity field in Figure 3. Notice that the inclusion of magnetohv- 
drodynamic effects, in contrast to the kinematic study of Low (1981). causes upward 
mass motion in addition to the up-lifting of the magnetic field lines because the plasma 
has to move with the field lines under the conditions of infinite conductivitv as 
manifested by the upward component of Lorentz force. Note, however, that some of 
the uplifted plasma (in the region displaced from the origin) slows down under the action 
of gravity, reverses direction, and falls back to the surface. Most of the motion, however, 
is upward. These upward mass motions are also found by Mikic. Barnes, and Schnack 
(1988) and Biskamp and Welter (1989). However, these workers did not include 
compressibility, pressure gradient, and gravitation as noted above. The present studv. 
which does so e.xplicitly, demonstrates a different evolution in the later stages. 

This induced upward motion can be explained via our governing equations. When 
we introduce the shear motion (r.), an axial field component. will be induced throuah 
the induction equation. The additional magnetic field will cause an additional maanetic 
pressure gradient in the momentum equation. This additional pressure gradient induces 
both the horizontal (i\.) and upward (i\.) motions as shown in Figure 3. Subsequenth. 
the mass motion interacts with both the magnetic field and gravity. Closer to the surface, 
the combined effect is dominated by gravity, and the result is the cluster of magnetic 
field lines in which a current sheet is formed as shown in Figures 2(g) and 2(h) at nearlv 
twice the Alfven time. 

Figure 4 shows the plasma properties (i.e., density temperature, and pressure en- 
hancement in terms of percentage change from the initial values at each level) at the end 
of this simulation (/ = 3600 s; more than 2 tj^). These properties are shown at various 
heights ( >’,, >'2, V4, >^6, and y,o, as shown in Figure 1(b)) as a function of horizontal 
distance. These results also help to explain the magnetic field line distribution. That is. 
the high density magnetic field region shown in Figures 2(g) and 2(h) within the 
mid-horizontal range (at the altitudes; ys, >"4) corresponds to the increase of plasma 
density by 20°; (i.e., dp/p ~ 0.2), temperature decrease of 20° „ (i.e..dr/T~ -0.2), and 
magnetic field strength (ABIBq) increase by a factor of 3. These properties are similar 
to those for a current sheet. With these properties in mind, let us now turn our attention 
to the plasma flow patterns as shown in Figure 3. The plasma flow rises initially above 
the zone of maximum shear velocity. At later times (say, from 1000 to 2000 s). the 
plasma flow moves toward the central region in a pattern reminiscent of a mushroom 
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600s 2400s 



1600s 3200s 



Fig. 3. Vectorial velocity field as a function of time during induced footpoint shearing motion for case (i): 

= 154; = 1700 s. 


cloud. In the later stages as shown in Figures 3(g) and 3(h), the significant plasma 
motion is again concentrated in the neighbourhood of the sheared region. This is also 
the region where the magnetic field lines have been clustered as seen in Figures 2(g) 
and 2(h). 
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Ap/p (p = 5.56 X 1 0" ^ g/cm^) 


AT/T (T =4500 k) 


A p/p (p =7.43x10^ dyne/cm^) 



Fig. 4. Changes ( relative to the initial local values ) of density {Ap = p - p^y{ y)), temperature {JT= T - T ). 
and pressure (dp ~ p ~ p(i(>’)) normalized by a proper value as shown at the end of the simulation (case (i); 

~ 154), t = 3600 s which is more than two Alfven time periods. The distributions are plotted alone the 
entire horizontal scale of the domain and at various levelst V| , va, and v’,,j as shown in Fieure 1(b). 
All the values are normalized by a reference quantity as indicated. 


2.2. Intermediate plasma beta (i.e., = 1.54) 

In this case, our simulation is performed with an initially modest magnetic field streneth 
(5o = 21.3 G) and with a shear velocity (wj of 15 km s ’ ‘ and - 46.5 km s ‘ The 
qualitative behaviour of the evolution of the vectorial fields (i.e., magnetic and velocity 
fields) and plasma parameters (i.e., density, temperature, and pressure) are similar to 
case (i). Therefore, we shall not repeat a full presentation. Nevertheless, there are some 
interesting features that appear in the evolutionary results of the magnetic and velocity 
fields as shown in Figure 5. 

The most pronounced result is the induced velocity distribution shown on the right 
side panels of Figure 5. The high velocity of the ascending movement in the central 
region is especially notable. As a result, the closed bipolar field tends to be opened up. 
We attribute this to the force created by the ascending movement of mass motion 
initiated by the shear prescribed at the lower boundary. The highest velocity attained 
by the mushroom cloud-like ascending mass motion is about 25 km s " ^ at r = 700 s 
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Po=1-54 
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Time = 100 s 


. , . ! t \ t . r~ 

I tl\ll / /»!.« 

|lt*l / 

til *1 \ % * • 

■ 0 t ^ ^ m 

^ ® ^ 

• ■ ,•* , • / ^ ^ m «*'N^ ,, 

I I . I 1 , I I . , , 
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Time = 450 s 




Fig. 5. Evolution of magnetic field lines and vectorial velocity fields at various times tor case (ii): fl, - 1 ^4. 
The characteristic Alfven time for this case is Tj, = 174 s. 


(i.e.. ~ 4 Alfven times) after introduction of the shear motion. The corresponding 
plasma parameters can be summarized as follows: the density decreases by about 50 
at the legs of the intermediate loops marked by the footpoints X 4 , X 5 . and .v^, as labeled 
in Figure 1(b). Again, the pinch effects discussed for case (i) occur and a current sheet 
is formed where the density increases by 25% ; the temperature decreases by 30°"o ; and 
the field strength increases by a factor of 2 . 
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2.3. Low PLASMA BETA (i.e.. = <J'06) 

In this case the initial magnetic field strength is increased to a more realistic value of 
106.3 G = 15 km s ' ' and = 232 km s ' ' without changing the other plasma 
parameters. The initial plasma beta is equal to 0.06 which is 250 times smaller than 
case (ii) and 2500 times smaller than case (i). .Again, the evolution of the magnetic field 
and velocity field exhibits patterns similar to those of cases (i) and (ii). Figure 6 shows 
the evolution of the magnetic field and the velocity vector field for this case. The 
maximum upward velocity is a factor of 4 higher than for case (ii) and a factor of 40 


Po=0.06 


Magnetic Field Line 


Velocity 


Time = 75 s 



%\\\\l\i////x 

x\vi\\\ 

ii \ / 

I . . I I I i I t I I . . 


(d) 


Time = 175 s 



Time = 225 s 



Fig. 6. Evolution of magnetic field lines and vectorial velocity fields at various times for case (lii): ~ 0.06. 

The characteristic Alfven time for this case is = 35 s. 
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higher than tor case (i). We note that the time required to reach the maximum velocity 
is much shorter than in the other two cases. 

In order to examine this phenomenon further, we plotted in Figure 7 the planar 
maximum absolute velocity (i.e.. (i\- + the neighbourhood of the apex of the 

arcade as a function of time for the three different cases. We choose to plot this 
parameter instead of the upward velocity. i',„ because the representative parameter 
[r; + r,-]‘ - is related to our analytical analysis that is discussed later (and in the 
Appendix). Actually, the numerical results show that the horizontal velocity. i\ . is only 
25% of the vertical velocity, l\.. First, we point out the change of scales that was required 
for the three cases (i), (ii). and (iii). Second, we direct attention to the common features : 
an approximately linear initial phase followed by a smooth transition to an explosive 
upward mass motion. The latter phenomenon is representative of the upward regions 
as discussed earlier. 


(a) 


Po- 154 



Time (s) 



Fig. 7. Maximum vectorial velocity that is representative of the upward vertical mass motion for cases ( 1 1 , 
(ii), and (iii). Note the change of scales. The representative Alfvcn times for the three cases = 154, 1.54. 
and 0.06, respectively) are = 1700 s, 174 s, and 35 s. 


It is interesting to relate these results to the magnetic field evolution. For example, 
we direct attention to Figures 2, 5, and 6 where, in the early stages of the evolution, the 
change of field lines is regular with a slowly ascending movement. This upward motion 
is also present in the force-free analyses of Low (1981) and Klimchuk and Sturrock 
(1989), and the numerical incompressible simulations of Mikic, Barnes, and Schnack 



364 


S. T. wi; ET AL. 


( 1988). and Biskamp and Welter ( 1989). However, the change of field lines in the present 
case becomes quite irregular in the later stages of the evolution. From Figures 2. 5, and 
6. we notice that the lower field lines are pinched together and the upper field lines tend 
to open up when the maximum planar velocity exceeds the Alfven speed. The Alfven 
speed for these three cases is 4.67 km s " 46.7 km s ' and 232 km s ', respectively. 

The maximum footpoint shear motion, v., is slow compared to the Alfven velocity in 
the latter two cases but fast compared with resistive diffusion in all three cases. Thus 
a sequence of essentially quasi-static, force-free states with frozen-in magnetic fields is 
found in the early stages, which ends when the magnitude of planar maximum velocity 
exceeds the Alfven speed, and the system becomes unstable. We claim that this is a 
shear-induced instability that could not be found in the earlier numerical simulations 
that omitted compressibility, pressure gradient, gravity, and the different treatment of 
boundary conditions. We shall return to this point later for further discussion utilizing 
analytical results. 

3. Further Interpretation of the Simulation Results 

From these simulation results, we have found that the buoyancy force leads to a 
mushroom cloud-like ascending movement that pushes the closed magnetic field up- 
ward. In order to understand this result further, we supplement our numerical simulation 
with an approximate analytical solution: 

3.1. Creation of mushroom cloud-like ascending motion 

From the numerical simulation of all three cases, we observe that the shear-induced 
mushroom cloud-like ascending movement can be ascribed to the out-of-plane com- 
ponent of the magnetic field. B.. This component gives an upward magnetic pressure 
gradient (i.e.. V(5f/8rt)) which causes the ascending movement of magnetic field and 
corresponding plasma flows. On the other hand, we notice that no 5. component is 
generated near the origin (.y = 0. r = 0) due to shear. This leads to a downward force, 
such that we observe the field lines being squeezed together to form a current sheet as 
shown in Figures 2. 5. and 6. This point can be illustrated further by using a linear 
approximation. The justification for the use of linear theory is seen from the numerical 
results that show that the initial stage of the shear-induced motion behaves regularly as 
shown in Figures 2, 3, 5, and 6. 

A closed form linearized solution for the induced field component B. is the following 
(for the derivation, see the Appendix); 

; — ■ = c 1 g cos(ax) cos[Lua (e cos(a.Y)) ' ] sin[(r + fo)Z.a>o ] ■ (3) 

This result expresses that the induced magnetic field B. rises from the lower boundary 
(i.e., y = 0) and spreads upward with a characteristic time scale LoJq, where L is defined 
by Equation (A. 8). It could be noticed from Equation (3) that B. decreases exponentially 
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with respect to the increase of v (height), because the term. cos[L< 2 .v \e ' cos(a.x)) ' ' ] 
in the central region, varies slowly with height. 

Finally, the coefficient c, corresponds to the shear velocity (n;.). The part of the total 
upward Lorentz force - ci cy{Bz ’Tf), that causes upward acceleration is 

independent of the sign of the coefficient c, (or w;.). 

3.2. She.xr-induced instability 

From the simulation results shown in Figure 7, we found earlier that instability sets in 
when the absolute maximum planar velocity exceeds the Alfven speed. In order to 
substantiate this claim, we performed a linearized analysis in which an approximate 
linearized solution for the planar velocities (u, i) was constructed as shown in the 
■Appendix (Equation (A. 13)). These velocities are as follows; 

u, = b' sin(2a.v) . 

r, = b' e ""'[1 -f- cos2(aa:)] . 

The electric current along the r-axis can be estimated, to the first order, as 

47t ^ _ cB, _ cB, 
c ' ex cy 

which means that the Lorentz force c “ - J^B.) leads to ascending flow, because 

it has been shown in the Appendix that b' is always positive and has an exponential 
growth rate as shown in Equation (A. 16). We have identified this phenomenon as the 
shear-induced instability since the numerical simulation results shown in Figure 7 are 
consistent with the analytical analysis. It is further noted from numerical results that 
the term - f “ 'J^B. is always upward. 

The results for the evolution of the magnetic field configuration shown in Figures 2. 
5, and 6 show clearly the two-stage evolution that we discussed earlier. The first stage 
of the evolution can be described by the linearized solution given in Equation (4). The 
seeond stage of the evolution involves the pinching together of field lines in the region 
where the shear motion was applied. If the three factors noted earlier (compressibility, 
pressure gradients, and gravity) had been absent, we believe that our results would have 
been similar to those of Mikic, Barnes, and Schnack (1988). Our current sheet, however, 
developed horizontally, whereas their current sheet was vertical. We explain this 
phenomenon by examining the distribution of upward component of the Lorentz force 
(i.e., c~ \J:-B^ - J^B.)). To illustrate this viewpoint, we use the results for = 0.06 
because this case best resembles the real physical conditions in active regions. The 
results are plotted in Figure 8. The left-most panels show the horizontal distribution of 
the vertical component of the Lorentz force at different heights from y^ to v,o (as show n 
in Figure 1(b)) at 25 s after the introduction of the shear motion at the lower boundary. 
A noted earlier, the Alfven time for this case is ~ 35 s. This result clearly indicates the 
first stage of the evolution due to the introduction of shear. All the forces are in the 
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t = 25s 


t = 100 s 


t = 213 s 


2.86 X 10 ^ dyn/cm ^ 



1.0G X 10'® dyn/cm 2 
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3.08x10 ' dyn/cm^ 
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Fig. 8. The total r-componcnt of the Loreniz force per unit area at / = 25 s. 100 s. and 213 s and at various 

levels in the solar atmosphere ( v = etc.). The representative .Alfven time for case ( lii) is 35 ^ md 

= 0.06. At r = 100 s (about 3r^) during the nonlinear stage of evolution, the Lorentz forces at the 
intermediate heights have a combination of upward and downward directions that causes magnetic held 
line pinching i see text ). This pinch effect is more pronounced at / = 2 1 3 s ( about ~ ) at lower altitudes. The 
horizontal axis represents the distance a'i . . . .v.. as shown in Figure lib) also shown for Figures 2-6. 


upward direction which means chat all field lines are lifted up in an orderly fashion. The 
magnitude of these forces is of the order of 3 x 10“*^ dyne cm “ The middle panels 
show the resultant upward component of the Lorentz force at r = 100 s which is about 
three Alfven periods. These results are reflected in the nonlinear nature of the evolution 
in which the Lorentz forces have both upward and downward direction at the inter- 
mediate altitudes. 

This bi-directional nature of the Lorentz forces causes the field lines to be pinched 
together in the lower regions as shown, for example, in Figure 6 for = 0.06. This 
particular feature is most pronounced in the results shown in the right-most panels 
which show the vertical component of Lorentz force at / = 213 s; this is about seven 
Alfven periods after the introduction of the shear. We note that the vertical component 
of this Lorentz force decreases at high levels, but, in lower levels (i.e., y, and y^), two 
very strong oppositely-directed vertical components of Lorentz force ( ~ 3 x 
10 “ ^ dyne cm ” *) appear. The force at y, is upward and the force at .vs is downward. 
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These iwo forces cause the held lines to be pinched together as shown in Figure o(c). 
Further discussion of this point will be included in the next section as part of a general 
scenario for shearing motions of magnetic arches or bipolar regions. 

4. Scenario 

From these simulation results, supported by the linearized analytical solution, a physical 
scenario is proposed for the formation of an 'Arch Filament System (AFS)' and its 
eruption as part of a more general scenario for 'Coronal Mass Ejections (CMEs) . A 
schematic representation of this scenario is presented in Figure 9. After introduction ot 
shear motion at a bi-polar region, all of the field lines will first be lifted up in an orderly 
fashion due to the shear-induced upward Lorentz force before the absolute maximum 
upward velocity reaches the local Alfven speed; this is the linear stage of the evolution. 
When this upward velocity is in the neighbourhood of the local Alfven speed, the lower 
parts of the magnetic field lines are pinched together, and an arch filament system is 
formed. At the same time, the upper part of the magnetic field lines is pushed upward, 
and a certain amount of mass is carried upward. This upward mass motion is shown 
in Figure 10 in terms of contours of dp and Ap that move upward at all but the lou esi 
gravitationally-bound heights. 



Fig. 9. Scenario for the formation of an arch-filament system (AFS) and upper level movement outward 
in the initial stage of a coronal mass ejection (CME) as a result of shear-induced instability. 
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Finally, when ihis absolute upward velocity exceeds the .Alfven speed, the shear- 
induced instability sets in as shown by the numerical results of Figure 7 and the 
analytical solution in the Appendix (Equation (A. 16)). In the following we compare this 
scenario with the available observations. 

Arch rilament systems and coronal mass ejections have been m\ esiigated b\ man> 
authors (Bruzek. 1967. 1968. 1969: Bumba and Howard. 1965; Martres ei uL. 1966; 
Harrison. 1986). These authors have noted that arch filament systems (.AFS) al\\a\s 
connect areas of opposite polarities and cross the neutral line in the longitudinal 
magnetic field. Bruzek (1969) has pointed out that the occurrence of AFS is associated 
with evolution of young bipolar spot groups. As for the motion of AFS. its characteristic 
feature is its expansion in height with an ascending velocity of 16-25 km s ' ' with 
footpoints rooted in the two opposite spot regions. This behaviour is quite similar to 
the early stage of the simulated magnetic field line evolution and mass motion shown 
in Figures 2. 3. 5, and 6 w here the apex of the magnetic loops is rising but their legs ha\ e 
little lateral movement. It was further noted that the AFS has both descending and 
ascending motions in loops. Bruzek (1968) attributed this phenomenon to the mass 
injection at one leg and its return to the chromosphere via another leg that has opposite 
polarity. On the other hand, shearing motion, if it has a line-of-sight component, would 
always lead to a blue shift in one leg and red shift in the other. Therefore, observations 
of flows in filaments arc not evidence of shearing. However, such evidence is not needed 
since the relative motion of bipolar spots is both necessary and sufficient evidence of 
shearing. Nevertheless, this concept of descending and ascending motion is based on 
Doppler shift measurements which can easily, at least partially, be recognized as 
complementary evidence of horizontal shear motion that occurs on both sides of the 
neutral line. This statement considers the fact that the spot group area is often not strictly 
perpendicular to the line of sight of the observer; thus the Doppler shift velocit\ must 
have an appreciable horizontal component (Harvey and Harvey. 1976). 

On the basis of our numerical simulations, the analytial solution and obser\cd 
characteristics, a physical model for the formation of AFS and subsequent CME can 
be constructed as follows. First, a young bipolar sunspot group emerges from the 
sub-photosphere. As it rises, its area increases and the neutral line dividing the opposite 
polarities gets longer and longer. Then a portion of the field can be reasonably regarded 
as a two-dimensional bipolar field (as is used in our mathematical model). In the 
meantime, the opposite polarity areas rotate with respect to each other. Associated with 
this rotation are horizontal shear motions that appear on both sides of the neutral line 
(thereby justifying our construction of the shearing velocity used herein). The Lorentz 
force generated by this process (see. for example. Figure 8) pushes the magnetic loops 
upward during an initial stage. .At the later times, the magnetic field becomes distorted, 
nonlinear MHD effects force field lines to pile-up and. then, the pinch phenomenon 
ensues. Such pinched magnetic flux tubes could be identified as arch filaments which 
are visible as a set of dark loops. The simulation has shown that in this region the 
plasma has high density and low temperature. From the analytical solution, we notice 
that the growth time ( F^a) " ' of the shearing instability is about 30 min which is a typical 
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Fig. i 1. Behaviour of ;/' and See Appendix (Equation (A. 12)). 


average life time of AFS. Thus, this simulation model may be appropriate to describe 
the formation of AFS and the eruption which leads to some CMEs. 

5. Concluding Remarks 

We have used a time-dependent, nonplanar MHD model for a bipolar maroetic resion 
that was subjected to shearing motion at its foot points. The characteristic plasma beta 
was varied over a wide range - from 154 to a more realistic value of 0.06. Common 
features were identified for all cases with the differences primarilv occurring in the timing 
of the events vis-a-vis the characteristic Alfven times. An essentially linear, early phase 
of upward mass motion was followed until the Alfven speed w'as reached, and a 
shear-induced instability is initiated. This nonlinear instability may be the basic 
mechanism for arch filament formation and subsequent coronal mass ejections. 

In our opinion, the early evolution in our simulation is in accord with quasi-staiic 
evolution of magnetic arcades demonstrated by Klimchuk and Sturrock (1989). In their 
work, a very low beta plasma was assumed, and therefore the magnetic field is unaffected 
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by pressure and graviiaiional forces. Our simulauons are also in accord with the 
dynamic evolution of magnetic arcades demonstrated by the numerical simulations of 
Mikic. Barnes, and Schnack ( 1988) and Biskamp and Welter ( 1989) in both the eariy 
and imennediaie stages of this evolution despite their neglect of compressibility, pressure 
sradient. and gravity. We did not rind the reconnection and formation of an ejected 
plasmoid. as Mikic. Barnes, and Schnack (1988) did. since we assumed electrical 


resistivity and viscosity to be zero. During the laie stages of the evolutionary develop- 
ment. when the plasma velocities surpassed the Alfven speed, our numerical simulations 
demonstrate nonlinear instability and catastrophic upward motion at high altitudes. 

As a final remark, it can be shown that these numerical results are valid over a 
wide range of parameters according to the scaling rule for dynamic similitude. For 
example, the present numerical results, computed on the basis of T,, = 10^ K and 
= 1.67 X 10~‘“gcm ■ ^. can be scaled to initial conditions of T, = 10*’ K ^nd 
p, = 1.67 X 10"''^gcm^ ’ by introducing a set of scaling parameters; ri = x ft- 


Z. I — t'l — X *■ 


r, = aT,,. p, = A- 'Po. P\ = Av and 5, = which leave the 


oovernina equations invariant for a given plasma beta. In a recent study ot similitude 


theory. Wu ei al. (1988) have shown that the present results also apply to the physical 
condition represented by these different initial conditions. 


As another example of the use of dynamic similitude, we may pose the following 
question: if the footpoints are moved slowly enough that the evolution is quasi-static, 
would the magnetic field closely approximate the static equilibrium states? Although, 
we suggested above (as did Mikic. Barnes, and Schnack, 1988, and Biskamp and 
Welter, 1989) that the answer is ‘yes', the reader is reminded of the values of the shearing 
velocity v, used in the present studies (e.g., 15 km s ^ maximum, for /5(, = 0.06) and 
in the above-mentioned work (30 km s ” assumed by Mikic. Barnes, and Schnack. 
1988. for ^ ~ 0.03). Although these maximum footpoint shearing velocities are much less 
than the Alfven speed, they are a factor of about 10 larger than observed photospheric 
velocities. 

In summary, we consider the results given here to be representative ot a realistic 
dynamical evolution of the posed physical problem of sheared magnetic arches and their 
evolution into arch filament eruption and coronal mass ejections. 

Finally, we remark on the relevance of our results to the observations of some CMEs 
as reported by Harrison (1986). The major point of his work is that a small X-ray burst 
is often found at the very onset of a CME, often followed by a large X-ray flare later 
on during the CME. In the present work, the formation of the current sheet coincides 
with the rapid increase in the velocity of the upper portion of the field lines. One could 
interpret the latter, as already discussed, as the onset of CME, while the current sheet 
formation could lead to a burst of energy dissipation (not shown here) which would be 
visible as a small X-ray burst. The simultaneity of these two events is consistent with 
the observations of Harrison (1986). This could be another indication that these numeri- 
cal results indeed represent a basic mechanism for the initiation of CMEs. 
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-Appendix 

To obtain an asymptotic solution for the relationship between the footpoint shearma 
velocity, ir. and B. in the first stage t linear stage) of evolution during which p. /;, T. . 
5,- vary slightly, we write 

P = Po+^>i, P=Pa^Pi. T=To+T,, 5, = . . 

( A. 1 ) 

5, = + 5,1 . B. = 5., , i\. = t\| , r, = r, , . r_ = r,, . 

where subscript 0 a nd 1 indicate the zero-order and first-order quantities. .And. i\ , 
r,.|i, ,f.,| Bni\ = .Alfven speed. B.y ' Bq. Inserting (.A.l) into Equations 
(2.4) and (2.7) formerly given by \Vu. Hu. and Nakagawa (1983) and leaving out the 
higher-order quantities, vve obtain the linearized equations 

^ g(g.i/y AttPo) ^ g,o 4^) _ h 5,0 5,, 

ct V 4ttPo cx ^ 4npr, cv 2 ^ 47rpo v 

c(5,,/^ 4;iPo) 5,0 ct., ^ 5 ,q cr,, 

ct y/47rpo c'.v y/4^ cr ’ 

where p„ = p e \ h = g/RT^. To solve Equation (A. 2), we construct the auxiliary 
equations 


ct 


cx 


cv 


c{B*!^Anp^ „ .cr: 

= 5,o(4;rpo) ' - ^ + 5,o(47rpo) “ ' ~ 

df t7.V CV 


(A.3) 


Substitutine 


\ s/4npoJ 

F* = . F~ = 


vf - 


y/47rpo/ 


2 


2 



Equations ( A. 3) reduce to 


cF^ 

ct 

cF~ 

ci 


5o(47rpj 


12 {a - h 2)y 


I 2 ^ - h l)y 


cF* cF^ 

cosax sma.t 


cx 


cv 


cF cF' 

-cosax + sina.x — 


C.Y 


Cl- 


(A. 4) 


Since solving Equations (A.4) is equivalent to solving their corresponding ordinary 
differential equations (Courani and Hilbert, 1962), it is easy to write down the solutions 
as follows; 


= (p(e - «'• cosa.Y, rcuo + fi^x) cosa.v) ' ^ ^ ^ 

f- = i//(e“"' cos<3.v. rcu„ - /(a,x) fe'"'' cos«.v)" ' . 

where 

\ 

cuo = aB„(4np^.)-' \ = j (cos .x' d.x' . 

0 

Considering the boundary value of c. (the nature of shearing) and using Equation i A.5) 
we can find the following solutions: 


r.* = c ^ cosax cos{Ll) s'miLrj) , 

( I = c, e~“'’ cosax sin(Z.iT) cos{Li]) , 

\v'47rpo/ 

where 

: = (r + r„)<^o • n = fi^x) cosax)~ ' ' . 

r„, I, and c, are integration constants. Back to solving Equations (A. 2) suppose i 
(5./^ 47rpo) satisfy the equalities (A. 6) except that L. c, are now not constants but 
functions of .x. r. Thus 


f., = c,(.x, >■) e cos(ajx) cos(L(.x, >-)l) sin(L(.x, r)t?) , 

[ f~‘ - ) = v) e""*’ cos(ax) sin(L(.x, y)0 cos(L(.x, y)r]) . 

\j4npoJ 

Inserting (A. 7) into (A. 2), c, and L can be determined uniquely by solving two ordinary 
differential equations. First, L satisfies the equation 


cosa.x sin ax — = ^(.x, y, L) , 

cx cy 

Q(x, y, L)= - (6/4) sin(ax) sin(2L 0 sin(2Z. rf) x 
X [( sin(2Z.j;) - sin(2LC)] ~ ' , 


(A.8) 
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with boundary condition L = L(x). After L has been found. (Inc, ) can be obtainea 
in the same manner using the following equation; 


cosifl.v) 


cflnc, ) 
cx 


sin(n.v) 


c(lnc.) 

— : = [y Ig(^0 - 7ctg(I?;)] Q{x. V. L) . 

cv 


I A. 9 1 


In fact, we only apply (A. 7) to explain the physical nature in the lower shearina reaion 
where ~ 0.8p^, therefore L and c, can roughly be regarded as constants. 

It is difficult to find an asymptotic solution for i\ and c,.. Let us consider case (iii) 
of strong magnetic field, in which the inenial force and - 7p and pg can safely be 
ignored. Inserting (A.l) into (2.2) and (2.3) of Wu. Hu. and Nakagawa (1983). the 
linearized equations are given as follows: 


cv 


Po 


vl 


B ' 




0 \ — I 

\ ov cx / 4 t: dx 


Po 


a 


47T V CV CX / 4n cv 


where the terms 


(A. 10) 


-L« 

4n cy 

which are second-order quantities, must be kept in view of actual mathematical manipu- 
lation. From (,A.7) the partial Lorentz force can be written as 


I r, cB., 
— 5., , 

47t ■ dx 


cB^ 

“ ~ '^ri = (cfn/2) (i/' - 1 - > 7 , ) c “ sin(2a.v) sin"(Z. 2) . 

cx 

(A.l 1) 
cB. 

“(47tPo )~ '^-1 — = {cfa/2)r]' + cos(2n.\)] sin-(Z,2), 

cy 

where ?; and i?, are slow-varying functions of .x. y. The representations for i; . p, are 
very complicated in the case with gravity, but we only deal with the lower central part 
of the domain where p~ const. Thus, the gravitational effects could be ignored in 
Equation (A. 2). then leading to the solution. n.xe“'' (cosa.x)' 5., ~ 5*. There- 
fore. r\' and p, asymptotically approach the case with no gravity. In such case p' and 
p, take simple forms as 


p' = {cos ny + n cos /r sin /j , 

p, = L cos n sin J7 (sin ax) ~ ‘ , 
n = L e‘^^'ax (cos a.x) ~ ' . 


(A. 12) 
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s. r. \vi et \l. 


Figure 1 1 shows the behaviour of t]’ and 17, . Note that if 1 is less than U.5. then 
0 < iji <1 17 — 1- Thereiore we will pay no attention to the difference between 17 and 
n + fh within the range ax: < (A.ll) reminds us of analogy between shearing 
velocity and force, so we suppose velocity having a mushroom-like form as 

I = sina.v . r,., = <3> - [ 1 + cos(2£/.v>] . (A. 13) 

where o' is a function oft. ,v. r (but weakly depends on ,v. r) being determined later. 
Inserting (.A. 13) into the linearized equations of (2.5) and (2.6). of Wu ei al. (1983) the 
time variation of current J.^ c can be found as 


cfcB^, cB .\ 

T ( ^ ) = I6a-5„d' ^ cos(fl.v) . (A 14) 

Cl V c.\ cr / 

In deriving Equation (A. 14) the weak dependence of d' on .y. r has been used. Differen- 
tiating (A. 10) with respect to i and inserting (A. 14) and (A.ll) into it and then letting 
it go to limitation when r goes to zero, we obtain one equation 

C‘~ 

— 0 ,_„ = 8r-;;a-(3 „ + {tTa/2)//' sin[2Z.aj„(r -t- ^,,)] (A. 15) 


to determine b' uniquely (here i-; = B^/4np,,). Noticing b'. //' only weakly depend on 
X, V, Equation (A. 15) can be regarded as an ordinary differential equation and. therefore, 
can be easily integrated with respect to t. Giving the initial condition: b' ! = 0. 

db /dr = 0 when t = 0. we obtain an asymptotic solution as 

^ Lv = o = [(2£ + ^)/2] exp(^/8 i ,^at) + [(a - (i)/2] exp( - ^ 8 v ^^ai) - 

sin[ 2 Z.o»o(/ + 6 ))] 

^ • r-,, 1 • (A. 16) 

sin[2Lcuor„] 

with 


x = 


La>r,ac^t] „ . Lr\'c}:i\, 

sin(2L£Uoro)- ' ^ 


8(L2cUo- + 2via-) 


8 y2 (Z." + 2) 


> 0 . 


R- li=0 , L~t]'c}iV^ 

P - Z — 7r cos(2LcUo^o) > 0 . 

16(L" + 2) 


8 \/2 -1- 2vXa~) 


Generally, we can find an approximate solution for the average b'. the representation 
of which is the same as (A. 16) except for the substitutions b' , = v = o. W , 

v^, where 


>■2 


b' = 


b' dv/vj , 


*1 = 


tj' dy/y- 


1x2 ’ 


vl = 



dy/y 2 
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From (A. 16) it can be seen that <3 will grow exponentially, and that the shearma veiocit\ 
c I acts like a 'seed'. If there is no 'seed', the mushroom flow velocities ( i\ . i\. j will never 
arise. The grow th rate is independent of c, but depends on the Alfven speed 
i\ = B,y^ Therefore, shear motion can induce linear MHD-instabilitv. However, 
this instabilit\ soon attains saturation, and the flow becomes quasi-steady and increases 
gradually until the velocities U\. ;\.) exceed i\. 
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